Lyapunov spectra of billiards with cylindrical scatterers: comparison with 

many-particle systems 
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The dynamics of a system consisting of many spherical hard particles can be described as a single 
point particle moving in a high-dimensional space with fixed hypercylindrical scatterers with specific 
orientations and positions. In this paper, the similarities in the Lyapunov exponents are investigated 
between systems of many particles and high-dimensional billiards with cylindrical scatterers which 
have isotropically distributed orientations and homogeneously distributed positions. The dynamics 
of the isotropic billiard are calculated using a Monte-Carlo simulation, and a reorthogonalization 
process is used to find the Lyapunov exponents. The results are compared to numerical results for 
systems of many hard particles as well as the analytical results for the high-dimensional Lorentz gas. 
The smallest three-quarters of the positive exponents behave more like the exponents of hard-disk 
systems than the exponents of the Lorentz gas. This similarity shows that the hard-disk systems may 
be approximated by a spatially homogeneous and isotropic system of scatterers for a calculation of 
the smaller Lyapunov exponents, apart from the exponent associated with localization. The method 
of the partial stretching factor is used to calculate these exponents analytically, with results that 
compare well with simulation results of hard disks and hard spheres. 

PACS numbers: 05.45.Jn 



I. INTRODUCTION 

Chaotic properties of systems with many degrees of 
freedom, such as moving hard spheres or disks, have been 
studied frequently. Extensive simulation work has been 
carried out on their Lyapunov spectrum 0, [E @]j and 
for low densities analytical calculations have been per- 
formed for the largest Lyapunov exponent 0, IE IE j the 
Kolmogorov-Sinai entropy |3, l9l 1 1 fj . and for the smallest 
positive Lyapunov exponents |lllll2j | . Many studies have 
also been done on the chaotic properties of billiards, sys- 
tems consisting of a point particle moving amongst fixed 
scatterers, notably the Lorentz gas [l3L ITU ITU lla| . 

In this paper, the effects of the shape and orientation of 
the scatterers in billiards are investigated, by numerical 
simulation and by analytical calculation. A comparison 
is drawn between systems of many freely moving hard 
disks or spheres and high-dimensional billiards with ran- 
domly oriented cylindrical scatterers. The hard-sphere 
system can be described as a single point particle moving 
in a high-dimensional space with fixed (hyper)cylindrical 
scatterers with specific positions and orientations. From 
the viewpoint of dynamical systems theory the Lorentz 
gas, where the scatterers are (hyper)spheres, is very sim- 
ilar to hard-sphere systems, as was noted already many 
years ago by Sinai |l7|. In a previous paper, we calculated 
the full spectrum of Lyapunov exponents of the high- 
dimensional dilute random Lorentz gas |13| . This spec- 
trum shows similarities to the spectrum of many hard 
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particles, but there also are differences. Some of these 
are due to the shape of the scatterers. Others result from 
their positions or their relative orientations. In this pa- 
per, these effects are further investigated by considering 
cylindrical scatterers. 

Section [n] serves as an introduction to Lyapunov ex- 
ponents. In Sec. IIIII the relevant dynamics of hard 
spheres are explained, and the correspondence to the 
hard-sphere system with cylindrical scatterers is dis- 
cussed in detail. In Sec. IIVI simulations are presented 
for systems consisting of homogeneously distributed, (hy- 
per) cylindrical scatterers with isotropically distributed 
orientations. These systems turn out to have spectra 
which closely resemble the spectra of hard disk systems. 
Section^contains a calculation of the smaller Lyapunov 
exponents of hard-disk systems. It relies on the recently 
introduced concept of partial stretching factors in 
combination with a carefully chosen isotropic approxima- 
tion. Finally, in Sec. IVII the results of both calculations 
are discussed further and a comparison is made. 



II. LYAPUNOV EXPONENTS 

Consider a system with a phase space T. In a d- 
dimensional system with N particles, the phase space 
consists of the positions and momenta of all particles and 
has 2dN dimensions. The system evolves from the initial 
state 7 at time t = 0, according to the path 7(70, t). 
If the initial conditions are perturbed infinitesimally, by 
J70, the system evolves along an infinitesimally different 
path 7(70, t) + 57(70, t), specified by 



(57(7o,t) = M 70 (i) ■ J 7o 



(1) 



2 



in which the matrix M 7o (t) is defined by 



,(*) 



(2) 



The Lyapunov exponents are the possible average rates 
of growth of such perturbations, i.e., 



Xi = lim — log i 

t— >oo t 



(3) 



where Hi(t) is the «-th eigenvalue of M 7o (i). The space 
of all perturbations of a point in phase space is referred 
to as the tangent space, denoted by ST. If the system 
is ergodic, the Lyapunov exponents are the same for al- 
most all initial conditions. The exponents are ordered 
according to size, with Ai being the largest and \2cLN the 
smallest, as is the convention. 

For Hamiltonian systems, such as hard spheres with 
only hard-core interaction, the dynamics are invariant 
under time reversal. For such a system, the attractor 
is invariant under time reversal, and so is the spectrum 
of Lyapunov exponents. Each tangent-space eigenvector 
which grows exponentially in forward time decreases ex- 
ponentially under time reversal. It is mapped onto an 
eigenvector with a corresponding Lyapunov exponent of 
equal size, but opposite sign. This is known as the conju- 
gate pairing rule. In systems such as the ones described in 
this paper, therefore, one only needs to calculate, numer- 
ically or analytically, the positive Lyapunov exponents. 
In systems which are reversible, but for which the attrac- 
tor is not invariant under time reversal, the conditions for 
and the form of the conjugate pairing rule are somewhat 
different 0. 



III. DYNAMICS OF HARD SPHERES 
A. In phase space and tangent space 

Vectors in tangent space which do not grow or shrink 
exponentially are generated by the symmetry operations 
associated with the symmetries of the dynamics of the 
system. They are eigenvectors with zero Lyapunov ex- 
ponents and are referred to as the zero modes. For a 
system of hard spheres under periodic boundary condi- 
tions, these symmetry operations arc uniform transla- 
tions, Galilei transformations, time translations, and ve- 
locity scaling. 

In order to calculate the remaining Lyapunov expo- 
nents, one must first derive the dynamics of the system 
in tangent space from the dynamics in phase space. Be- 
low, the position and velocity of of a specific particle, 
indexed by i, will be denoted by r, and Vj. 

The evolution in phase space consists of an alternating 
sequence of free flights and collisions. During free flights 
the particles do not interact and the positions grow lin- 
early with the velocities. At a collision, momentum is 




FIG. 1: Geometry of a collision of two particles of diameter 
a in relative coordinates. The collision normal a is the unit 
vector pointing from the center of one particle to the center 
of the other. The circle delineates the collision surface, the 
distance of closest approach. 



exchanged between the colliding particles along the col- 
lision normal, a = (r, — Vj)/a, as shown in Fig. ^ The 
other particles do not interact. 

From these dynamics, the tangent-space dynamics can 
be derived. During free flight there is no interaction be- 
tween the particles and the components of the tangent- 
space vector transform according to 



6v' 



1 (t-to)l 
1 



<Sv ? ; 



(4) 



in which 1 is the d x d identity matrix and the primes 
indicate the vectors after the free flight. 

At a collision between particles i and j, only the 
tangent-space vectors of the colliding particles are 
changed |£g. As shown in Fig. ^ an infinitesimal dif- 
ference in the positions of the particles leads to an in- 
finitesimal change in the collision normal and collision 
time. The v + 5v are exchanged along a + 5a. This 
leads to infinitesimal changes in both the positions and 
velocities right after the collision. For convenience we 
change over to the relative and center-of-mass coordi- 
nates, Srij = Sri — 5r , . <JR , , = (<5r; + 5vj)/2, Vy 
v. - Vj,Vij = (vj + Vj)/2,Svij = 5vi - 5vj, and 
SVij = (Svi + Svj)/2. The tangent-space coordinates 
are found to transform as 



(5r' = 8ra -25 -Si 



SB.L 



SR,. 



8V l3 = 8v i:i - 2S ■ 8v i:i - 2Q • Si 
tfV' = SVu , 



(5) 
(6) 
(7) 
(8) 



in which S and Q are dxd tangent-space collision matrices 
and primes are used to denote the coordinates after the 
collision. The matrix S may be written as 



S = 



(9) 



3 



where the notation ab denotes the standard tensor prod- 
uct of vectors a and b. Define 9 as the angle between v i:) - 
and <7. The unit vector orthogonal to Vy in the plane 
spanned by and a then reads 

(1 — V;,-V„) - a , 

9 |sin0| 1 ' 

The matrix Q may be written as 

Q = ^ (cos6(l -VijVij - pp) + -!— p'p) . (11) 
a \ cos ft J 

Note that Q transforms the vectors Sr^ which are or- 
thogonal to Vij into vectors that are orthogonal to Vy. 
The vector v,j is a right-zero eigenvector of Q, and Vy a 
left-zero eigenvector. Note that these are d-dimensional 
vectors, not 2(i-dimensional. 

B. Scatterer configurations 




FIG. 2: The two-dimensional representation of a system con- 
sisting of three particles in one dimension. A possible trajec- 
tory is shown. If the system is infinite, the particles will soon 
leave each other's vicinity and never collide again. The axes 
of the cylinders intersect in one point. If the system is subject 
to periodic boundary conditions with periodicity L, they will 
encounter each other again. The cylinders belonging to such 
collisions and a path are indicated with dotted lines. 



In the relative coordinates of the two particles, the col- 
lision surface is spherical. As is shown in Fig.^ the rela- 
tive coordinates transform as the coordinates of a point- 
particle colliding with a sphere of the same dimension. 
The center-of-mass coordinates do not transform. 

The dynamics of two colliding hard disks or spheres 
are equivalent to the dynamics of a point particle collid- 
ing with a fixed (hyper)cylindrical scatterer. A system of 
many hard particles thus is equivalent to a system con- 
sisting of a point particle bouncing between many cylin- 
drical scatterers. Each of the scatterers corresponds to 
the collision surface of two particles. 

The orientations of the cylinders can be readily found. 
If two cylinders belong to collisions of pairs of different 
particles, the finite directions of the two cylinders are 
orthogonal. Two specific cylinders may be collision sur- 
faces belonging to two pairs of hard spheres which involve 
a common particle i. In this case, the spherical compo- 
nents of the two cylinders are not orthogonal. Let the 
two other particles involved in the two collision surfaces 
be particles j and I. The (i-dimensional plane with which 
the intersection of the collision surface of particles i and 
j is a sphere then consists, if I ^ i,j, of the sets 



Si- 



{r|r. t = -Tj Ar; = 0} 



(12) 



For any d, the highest possible value for an inner product 
of unit vectors in <Sjj and Su is ^ ■ Therefore the angle 
between the two sets is 7r/3. 

The coordinates of the center-of-mass constitute d out 
of the dN dimensions of the original system. If the 
boundary conditions are periodic or if the system is infi- 
nite, there is only uniform motion in these directions. As 
a consequence, perturbations in these coordinates (uni- 
form translations of the entire system) remain constant. 

Perturbations of the velocity that correspond to Galilei 
transformations lead to a linear growth in the perturba- 
tions in the positions. This yields 2d linearly indepen- 
dent perturbations which do not grow exponentially, and 



therefore 2d zero Lyapunov exponents, associated with 
the position and momentum of the center of mass. There 
are two more zero Lyapunov exponents, corresponding to 
a translation in time and a rescaling of the velocity in dN 
dimensions. 

If the system has periodic boundary conditions, there 
are extra cylinders, corresponding to collisions after at 
least one of the particles has moved through the bound- 
aries of the periodic volume. After eliminating the center- 
of-mass coordinates by setting the center-of-mass posi- 
tion and momentum to zero, one can describe the sys- 
tem as a point particle with coordinate r in a d{N — 1)- 
dimensional space with periodic boundary conditions 
moving among N(N— l)/2 fixed (hyper)cylindrical scat- 
terers with d spherical dimensions. As energy is con- 
served, the particle still moves at velocity v, which is 
related to the inverse temperature (3 = l/{k-oT) and the 
particle mass m by 



[N — l)d 



(3m 



(13) 



Consider, for example, the simplest nontrivial case of 
three particles in one dimension. The space has two di- 
mensions and the fixed (hyper)cylinders have one spher- 
ical and one infinite dimension. The 2-dimensional rep- 
resentation is displayed in Fig. [5] 

The cylinders corresponding to the hard-sphere sys- 
tem are oriented in specific directions. To simplify cal- 
culations, it is possible to consider a homogeneous distri- 
bution of scatterers, with a distribution of orientations 
which is isotropic. With such a distribution and the ap- 
proximate distribution of the radius of curvature tensor 
derived in Ref. @, it may become possible to use the 
techniques developed in Ref. 01 to calculate the Lya- 
punov spectrum. 



4 



IV. SIMULATIONS OF THE SPECTRUM OF 
ISOTROPICALLY DISTRIBUTED CYLINDERS 

To investigate the effects of the shape of the scatter- 
ers and the distribution of orientations, simulations have 
been carried out for a high-dimensional system with ho- 
mogeneously distributed cylinders with an isotropic dis- 
tribution of orientations. Cylinders are considered with 
two spherical directions. This system can be compared 
to hard disks with the same collision frequency as well as 
to the Lorentz gas. 

A. Simulation method 

The dynamics in phase space have been calculated by 
means of a Monte Carlo method. Collision parameters 
were drawn from the relevant distributions, to generate 
a path in phase space. Each step consists of a short free 
flight of the point particle, followed by a collision with a 
cylindrical scatter. The free-flight times are drawn from 
an exponential distribution of times r, 

p(r) = v N exp(-%-r) , (14) 

where Dn is the average collision frequency of the point 
particle with the scatterers. We have 

v N = ^DN , (15) 

with v the average collision frequency of a single particle, 
at low densities given by 

_ _ 2^- 1 )/ 2 na rf ~ 1 

r(rf/2)v^ ' (1 ' 

The collision frequency of a specific particle depends 
on its velocity. The collision normal is drawn from an 
isotropic distribution in the d(N — l)-dimensional space 
such that a ■ v < 0. It is accepted with a probability 
equal to the size of the component along the velocity 
of the point particle, —a ■ v, which is proportional to 
the collision rate in two dimensions. The orientation of 
the scatterer is specified by the two vectors in its spheri- 
cal directions, the collision normal and one other vector, 
which is drawn from the isotropic distribution as well. 
This leads to an isotropic distribution of the orientations 
of the scatterers. For more details on Monte Carlo sim- 
ulations, see Ref. [is||. 

At each collision, the dynamics of the point particle 
are calculated, as well as the transformations of a num- 
bered set of tangent-space vectors. After every step, the 
tangent-space vectors are reorthonormalized. That is to 
say, the components of each of the vectors along vec- 
tors with higher indices are subtracted and then the vec- 
tors are normalized. The scaling factors are equal to the 
growth of each vector between the last two collisions. For 
long times, the growth of the i-th vector is dominated 



by the i-th Lyapunov exponent. The scaling factors are 
stored and, for each vector, their logarithms are summed. 
The i-th sum divided by the total elapsed time converges, 
for long times, to the i-th Lyapunov exponent. For more 
details on this method for calculating the Lyapunov ex- 
ponents, see Ref. |l9j . 



B. Discussion of the spectrum 

One may compare the spectra of a point particle col- 
liding with homogeneously and isotropically distributed 
cylinders to the spectra of systems of hard disks with 
the same collision frequency, energy, and dimensional- 
ity (d = 2), as well as to the spectrum of the high- 
dimensional Lorentz gas. Plots of the spectra of the 
isotropically distributed cylinders at various dimensional- 
ities and densities are shown in Fig. [31 The spectra for the 
corresponding hard-disk systems with periodic boundary 
conditions are shown in Fig. 01 

The Lyapunov spectrum of the high-dimensional 
Lorentz gas, calculated in Ref. i s much flatter than 
that of the hard-disk system. It becomes flatter with 
increasing number of dimensions. This difference in be- 
havior is due to the fact that, for the Lorentz gas, all 
coordinates not associated with zero modes are involved 
in every collision and grow with a factor of the order of 
the free-flight time between two collisions. For hard disks 
only the four (2d) phase-space coordinates of the two col- 
liding particles are involved in a collision. From Figs. EH 
and 01 one can see that the spectrum of a point parti- 
cle colliding among isotropically distributed cylinders is 
much more similar to the hard-disk spectrum. For large 
iV", its shape becomes independent of N. 

For the largest exponent it is known that the corre- 
sponding perturbation is carried by only a few particles 
|4|. The Lyapunov exponents are strongly affected by 
this, because collisions of particles other than these few 
do not contribute to the growth of the tangent-space vec- 
tors. There is a small probability of large growth, as op- 
posed to a large probability of small growth in the fully 
isotropic system. In the lower end of the spectrum, the 
tangent-space vectors for hard disks are carried by many 
particles |2d| |. In this regime the exponents behave simi- 
larly. From Figs. |3| and ^ it can be seen that they depend 
differently on both the density and the particle number. 

The tangent-space eigenvectors corresponding to the 
smaller Lyapunov exponents of hard-disk systems are 
carried by many particles, and for these the similarities 
to the cylinder system are much greater. The lower Lya- 
punov exponents behave similarly to the exponents of the 
hard-disk system (apart from the Goldstone modes) and 
are proportional to v. For the d = 2 cylinder systems, 
I find that the smallest positive Lyapunov exponent be- 
haves as 

A 2 (at-i)-2 = 0.37z/ . (17) 
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FIG. 3: The spectra of positive Lyapunov exponents for sys- 
tems of isotropically distributed (hyper)cylinders with the 
same collision frequency, energy, and dimensionality (d — 2) 
as hard-disk systems with densities 0.1a -2 , 0.01a -2 , 0.001a -2 , 
and 0.0001a -2 , and various particle numbers N. The inverse 
temperature j3 = 1. 



For large particle numbers and low densities, the small- 
est exponents of hard disks an d sp heres, apart from the 
Goldstone modes, are equal to |21j 



A 



d(JV-l)-2 



0.31 v 
0.39 v 



if d = 2 
if d = 3 



(18) 



There is a difference of less than 20% between the small- 
est exponents of the cylinders and the hard-disks, inde- 
pendent of density or particle number, provided the par- 
ticle number is sufficiently large and no Goldstone modes 
exist in the hard-disk system. 

The Goldstone modes in the hard-disk system (see, for 
example, Refs. jHHHlllllil) are due to localization. 
In the corresponding high-dimensional system of cylin- 
drical scatterers, they are related to the specific positions 
and orientations of the scatterers. As the cylinders in the 
simulation presented here cannot be associated with two 
specific particles, Goldstone modes are absent. 



FIG. 4: The spectra of positive Lyapunov exponents for sys- 
tems of many freely moving hard disks, with the same collision 
frequency, energy, and dimensionality (d — 2) as the systems 
for which the Lyapunov spectra are displayed in Fig. [3] 



V. ISOTROPIC-CYLINDER 
APPROXIMATIONS FOR THE HARD-DISK 
SYSTEM 



The simulation described in the previous section pro- 
duces a Lyapunov spectrum which is suggestively similar 
to the spectrum of hard disks. As the system is homoge- 
neous and isotropic, the techniques developed in Ref. [ijj 
can be applied to calculate the Lyapunov exponents of 
the isotropic cylinder system. Rather than working this 
out, one may return to the hard-disk system and apply 
some approximation inspired by the similarities between 
the two systems to the calculation of the Lyapunov ex- 
ponents of hard disks. In the lower end of the spectrum, 
the tangent-space vectors for hard disks are carried by 
many particles |2f| . and in this regime it is possible to 
use the isotropic approximation to calculate Lyapunov 
exponents. 

In this section, a calculation of the smallest exponents 
not due to localization effects is described which makes 
use of the stretching factor derived in Ref. [8j , but which 
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assumes an isotropic distribution of scatterers. This ap- 
proach will lead to behavior of the largest exponents that 
is different from the real hard-disk system. However, for 
the smaller Lyapunov exponents the approximation will 
be better. The Goldstone modes, which were discussed 
in, for example, Refs. 0, H2, L^J 123, l2J| > are also removed 
by the homogeneous distribution of scatterers, as their 
existence relies on the fact that in the hard-disk system 
only nearby particles collide. 

A. Partial stretching factors 

This subsection briefly discusses the method of the par- 
tial stretching factor [13j to calculate the full Lyapunov 
spectrum of an isotropic system. In standard terminol- 
ogy, the stretching factor is defined as the factor by which 
the expanding part of tangent space expands over a time 
t. It can be used to calculate the Ruelle pressure as well 
as the sum of the positive Lyapunov exponents, equaling 
the Kolmogorov-Sinai entropy in systems without escape 

nam. 

By analogy, the partial stretching factor As(r, v, t) of 
a p-dimensional subspace S of the 2cW-dimensional tan- 
gent phase space is defined as the factor by which the 
volume of an infinitesimal p-dimensional hypercube in 
this subspace has increased after a time t. Unless S is 
orthogonal to one of the p most unstable directions in 
tangent phase space, the partial stretching factor will, 
for very long times, be dominated by the p largest Lya- 
punov exponents. Explicitly, one has the identity 

p i 
VA,= lim -logA s (r,v,*) . (19) 

— ' t^oo t 

4=1 

The partial stretching factor just after collision N is the 
product of the partial stretching factors due to collisions 
1 through N. These depend on the relative orientations 
of v, a, and the image of S. 

For systems with only hard-core interaction, where the 
collision times are exactly defined, the partial stretching 
factor can be written as the product of partial stretching 
factors resulting from each of the different single colli- 
sions combined with the subsequent (or previous) free 
flights of the two particles involved. In this description, 
the effects of the free flights of the other particles have al- 
ready been accounted for at their most recent collisions. 
On the right-hand side of Eq. (|19|) . the logarithm may 
be replaced by the sum of logarithms of these stretching 
factors. The resulting expression may be interpreted as a 
time average, which in ergodic systems may be replaced 
by an ensemble average. If the system is isotropic, the 
collisions are uncorrelated, so that one has 

± A i= ^(logAf) , (20) 

4=1 

which is independent of the chosen initial subspace. Here, 
A<f ] is the single-collision partial stretching factor due to 



collision j of a p-dimensional subspace of tangent phase 
space. This is assumed to include the free flights of 
the particles after the collision, and not those before. 
To obtain the Kolmogorov-Sinai entropy and the Lya- 
punov exponents, one must calculate the distribution of 
single-collision partial stretching factors. For more de- 
tails on the derivation of Eq . (|20|l and the consequences 
of isotropy, see Refs. jlj, |25| . 

The growth of a d/V-dimensional volume element in ST 
can be monitored through its projection onto a subspace 
of ST with at least the same number of dimensions, as 
long as this projection space is not orthogonal to one of 
the dN leading eigenvectors of M. In the limit t — > oo, 
the logarithm of the determinant of the transformation of 
the projection yields the same Kolmogorov-Sinai entropy 
as the logarithm of the stretching factor of the actual 
volume clement. 

If (5r( m \ Sv^) are the eigenvectors belonging to the 
positive exponents, the eigenvectors which belong to 
their counterparts under conjugate pairing are equal to 
(Sr^ m \ — cSv 1 ™)). Therefore, eigenvectors which have no 
contributions along either Sr or Sv correspond to them- 
selves under conjugate pairing. Such eigenvectors must 
therefore have zero Lyapunov exponents. The spaces 
spanned by either Sr or Sv are not orthogonal to any 
eigenvectors which belong to nonzero Lyapunov expo- 
nents. In the system described here, a convenient choice 
for the projection space may therefore be either of these 
spaces. Often, Sv is used, because it does not change 
during free flights. However, in the calculation presented 
here it is necessary to use Sr instead, as I will show in 
Sec. EH] 



B. The inverse radius of curvature 

This subsection briefly summarizes the relevant results 
and intermediate results of Ref. 8], from which more 
details can be obtained. 

In order to calculate the transformation of the pro- 
jection of a tangent space vector onto a subspace of the 
tangent space, one needs information about the original, 
unprojected vector. The single-particle partial stretching 
factor depends on Sr as well as Sv before the collision. 
As dN dimensions will be projected out, Sr may be as- 
sumed to have a probability distribution which depends 
on Sv. The radius of curvature is defined as the tensor 
which transforms Sr into <5v. The inverse of the radius 
of curvature connects the two perturbations, 

Sr = ?W-Sv, (21) 

with f — 1/v the average free-flight time. The matrix 
W can be split up into d x d matrices between specific 
particles, Wy. If the two indices are equal, Wj may be 
used as shorthand. As particles collide and have free 
flights, Wjj changes. The volume element projected onto 
Sr or S~v before the collision is mapped to a projection of 
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a volume element after the collision. This map depends 
on the elements of W. 

Maps of W can be found from the dynamics. The 
distribution of elements of W may not change by a colli- 
sion. Together with the distribution of collision parame- 
ters this yields a complicated nonlinear integral equation 
for the joint distribution function of the elements of W. 
With p(W) the distribution of W before the collision and 
p'(W) the distribution after the collision, one may write 

p'(W) = J dW P {W)s(w'(W)~W^ . (22) 

In Ref. the equation is solved by an iterative 
method, starting from a fairly simple distribution and 
iterating the equation. Each extra step in the iteration 
results in a distribution function which more closely re- 
sembles the true solution. In order to improve the con- 
vergence, I introduce a parameter, the average trace ele- 
ment, in such a way that the average trace can be kept 



fixed over an iteration. In the first iteration, just before 
a collision, the are equal to their averages, and, due 
to the Stofizahlansatz, all \Nij are equal to zero. If the 
distribution of the angle between the relative velocities 
of two consecutive collisions is (nearly) isotropic, the two 
average diagonal elements are (approximately) equal. In 
this case, the matrix is 

Wy = wlS tJ , (23) 

where Sij is the Kronecker delta. The initial distribution 
used in the iteration process is a product of Dirac delta 
functions at the average value w for the diagonal elements 
and zero for the off-diagonal elements. 

From the dynamics, formulated in Eqs. |@J|-||SJ|, one 
finds W after the collision (W*) and after the collision 
and the free flight (W). In the basis consisting of 
and the d — 1 vectors orthogonal to it, the values of Wjy 
are changed according to 



{w) 
{\w)U-i 

\ 
-hiol d -x ) 



if k = l = iVk = l= j , 



if (k,l) = (i,j)v(k,l) = (j,i) , 



(24) 



wlS, 



k-l 



if k^i,jVl^i,j , 



(w + VT k ) 





(\w + VT k )\ d -l 



if k = I = i V k = I = j 



w, 



kl 




1 

2 

w 1 S k i 



-hwld-x 



if (k,l) = (i,j)V(k,l) = (j,i) , 



if k ^ i,jV I ^ i,j 



(25) 



where ld-i denotes the (<£— l)-dimensional identity ma- 
trix. Equations l|24|) and (|25|l imply a distribution for the 
elements of expressed in the basis belonging to the 
next collision, which consists of v£ - and the d — 1 vectors 
orthogonal to it, vLj_- 

After two iterations of Eq. (|22p. it is found that the 
average trace element is approximately equal to 



w 



2.929 if d 
1.947 if d 



2 , 

3 . 



(26) 



Together with Eq. (|25|l and the distribution of the basis- 
transformation matrix, which may be conveniently drawn 
from a molecular dynamics simulation of hard disks or 



spheres, this gives an approximate joint distribution func- 
tion for the elements of W,- . 



C. Projection 

For the Lorentz gas the isotropic distribution of scat- 
terers makes it possible to simplify the calculation [l3| . 
Because of the isotropy, the probability distribution 
of the stretching factor is independent of which p- 
dimensional subspace S is being stretched. Also, the 
low-density approximation is not affected by the prob- 
lems described in the previous subsection and in Ref. |j| . 

For systems composed of isotropically distributed 
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cylindrical scatterers or hard disks, however, the choice 
of the space onto which everything is projected affects 
the distribution of the partial stretching factors. Diag- 
onal elements of W grow linearly during free flights. In 
the calculations in hard-disk systems in Ref. 0, the el- 
ements of W are, in fact, calculated as weighted sum- 
mations over sequences of free flights. Notice that in 
the case of truly isotropically distributed cylinders the 
summations are somewhat different, and that the partial 
stretching factors are therefore different from those in the 
equivalent hard-disk system. 

The free flights of several collisions contribute to an 
element of W, as with each free flight Sr grows linearly 
with Sv. This can be seen, for example, in Eq. (|25Jl . 
Suppose one such term contains a long free-flight time. 
In the spatial tangent-space direction belonging to this 
particle, the stretching is large, and this affects the orien- 
tation of the stretching manifold. If the growth is mon- 
itored through a projection onto Sv, the stretching due 
to several free flights affects the partial stretching factor 
at one collision. At another collision of one of the two 
particles, the same long free flight will again affect the 
partial stretching factor. This introduces a correlation 
between the orientation of the projection of S onto 5v 
and the amount of stretching. The distribution of par- 
tial stretching factors then depends on the orientation of 
the stretched space S, even in the case of fully isotropi- 
cally distributed cylinders. In such cases, Eq. (|2(J[1 cannot 
be used. 

If an approximation not exhibiting this problem is to 
be made for the hard-disk system, the standard choice of 
projection on <5v is inadequate. Instead, as is done in the 
present paper, the dynamics must be projected onto <5r. 
In this representation, the correlations through the sums 
of collision times are removed, because the linear growth 
is incorporated in the partial stretching factor immedi- 
ately. However, the isotropic distribution of orientations 
still is an approximation, as there is correlation between 
free-flight times and other collision parameters from dif- 
ferent collisions through the particle velocity. 



D. The partial stretching factors for hard spheres 
at low densities 



ing factor of the projection onto Sr after a collision and 
subsequent free flights is, to leading order in the density, 
2vt + I (a cos 6) . Here, r + is the same as in Ref. @, that 
is, r + = (ji + Tj)/2. In d — 2 directions orthogonal to 
p and v,j, the partial stretching factor is 2vt + cos9/a. 
There are d + 1 coordinates involved in the collision on 
which Q does not work, the center-of-mass coordinates 
and the relative coordinates parallel to v.y. In these di- 
rections, the linear stretching due to the free flights must 
be calculated and incorporated into the partial stretch- 
ing factor. Meanwhile, for all other particles, Sr grows 
linearly as well. This can be accounted for later, at their 
next collision, without loss of generality. 

The distribution of the stretching is difficult to obtain 
if the distribution of the elements of W is complicated. 
In a simple calculation, the expressions for W' and W* in 
the approximation of Sec. IV Bl can be used [see Eqs. 
and (I25|) ]. In this case, W is approximated by W = 
lw. With this simple form of W' and W*, the partial 
stretching factors in the remaining d + 1 directions can 
be calculated. 

For the d — 1 directions of the center-of-mass coordi- 
nates orthogonal to Vy, the eigenvalues of W are w + 
Pt + , yielding a partial stretching factor of (w + vt + )/w. 
Similarly, for the remaining two directions, those paral- 
lel to Vij, one finds (w + uTi)/w and (u> + vTj)/w. In 
the directions belonging to particles not involved in the 
collision, nothing changes. In the coordinates parallel to 



• ^ = (l + f ) 



Sr, 



(27) 



In the relative coordinates, Q acts on vectors orthogonal 
to v'„ , 



(1 - Vyv^-) • Sr' i:j 
2vt. 



cos 0(1- vLvy -p'p) 



1 



cos 9 



p'p 



Sri 



+ £CZL_2) (1 _ w .fl^ 



(28) 



The stretching of the projection onto Sr of a p- 
dimensional subspace of the tangent space due to a colli- 
sion between particles i and j depends on this projection 
and on W at the collision. There are 2c? coordinates 
involved in the collision projected onto <5r, d for each 
particle. The tangent-space dynamics are described in 
Sec. IIIII In the relative coordinates, the collision trans- 
forms the tangent-space vectors as in an elastic collision 
with a G?-dimensional, fixed, spherical scatterer. Q works 
on d — 1 directions of the relative coordinates orthogonal 
to Vy . The action of Q on the relative coordinates is 
described by Eqs. Q and (|TT)) . 

In the direction of <5ry parallel to p, the partial stretch- 



In the center-of-mass coordinates orthogonal to vL, 



(l-v^)-^=(l 



(1 - yljVij) ■ SR 



2w 



-(1- VyVij) -5r 



(29) 



The eigenvalues of the transformation, the growth fac- 
tors, are denoted by gi, with / between 1 and 2d. They 
can be found from Eqs. I|27|l - (|29|l . In summary, the 
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following growth factors occur: 



91 



2vt + 


if 


1 = 1, 


LI LUo U 






2vt + cos 9 


if 


Kl<d-1, 


a 






VTi 
1 + — 


if 


l = d, 


w 




l + ^± 


if 


d < I < 2d - 1 


w 






1 + ^ 


if 


I = 2d . 









(30) 



From these growth factors, the partial stretching factor 
can be calculated for any subspace S of tangent space. 



E. Lower bound of the Kolmogorov-Sinai entropy 

In principle, the choice of projection does not affect 
the Kolmogorov-Sinai entropy, since its calculation does 
not involve partial stretching factors, only the stretch- 
ing factor. There is, therefore, no need for an isotropic 
approximation. From Eqs. I)26|l and l|30|l the stretching 
factor can be calculated. From Eq. I|20|) the Kolmogorov- 
Sinai entropy is found to satisfy 



h^s = — ( log 



' 2d ' 

i=i . 



(31) 



The approximation of W^ 1 by 1/VV affects the 
Kolmogorov-Sinai entropy. After numerical integration, 
using the estimate for w in Eq. I|26|l . this yields for the 
constant B in the expansion hits = AND[— log(na d ) + 
B + . . . ] the approximate values 



0.98 if d = 2 
0.13 if d= 3 



(32) 



The estimation obtained here is less accurate than the 
results of Ref . Q and the numerical values found in sim- 
ulations Hj. The latter are 1.366 and 0.29. This is 
related to the fact that the distribution of the stretching 
factor after only one iteration of the equation for the dis- 
tribution function is used. Note that a wider spread of 
the elements of W, or a lower value of w, which would 
result from more iterations, leads to a larger value for B. 



F. The smallest exponents 

Despite the approximate nature of the calculation of 
the partial stretching factors in this section, it is possi- 
ble to use the results for an estimation of the smallest 
Lyapunov exponents. 

Any [d(N — 1) — 2]-dimensional subspace £d(jv-i)-2 
of the [d(N — 1) — l]-dimensional subspace Smn— of 



[6ri] orthogonal to the zero modes can be characterized 
by a single vector, g. This is the vector orthogonal to 
Sm jv-i)-2 an d the zero modes. The approximation made 
in this section is that this vector, which is the eigenvector 
belonging to the smallest positive exponent, has signifi- 
cant com pon ents along the directions Sri and 5vi of many 
particles pCJ | , and is more or less isotropically distributed 
in phase space. Also, these components do not depend 
much on the velocities of the particles. This permits the 
isotropic approximation. 

The smallest exponent can be calculated from Eq. l|19fl . 
the partial stretching factor of S<j(jv-i)-2 a t a collision, 
and the stretching factor [l^. Following the derivation 
in Ref. for the hard-disk system, one obtains 



*d(N -l)-2 — \ l0 S A i _ l0 S A i 



(33) 



where i is the index of the collision, Aj is the stretching 
factor due to collision i, and A^ ci ' A '~ 1 ' ) ~ 2 ' > is the partial 
stretching factor of Sd(N-i)~2 due to collision i. 

The difference of the two logarithms in Eq. I|33fl can 
be expressed in terms of the components of g along the 
growing directions, denoted by sin^/, with I the index of 
the growth factor. If sin</>/ is small for all I, (sin 2 cf>i) — 
l/(dN). One finds 



2,/ 



Ad(AT-l)_ 2 ~ ( - log ( Y[ \l COs2 & + ~2 Sin2 ' 



. 1=1 



(34) 



If g has significant components along many particles, 
then sin <f>i is small, and the logarithm can be expanded 
about unity argument, to yield 



A 



d(JV-l)-2 



EM 



(2=1 



sin 2 4>i 



(35) 



If <f>i is isotropically distributed in d(N — 1) — 2 dimensions 
and N is large, the average of sin 2 (f>i can be approximated 
by l/(dN) + 0(l/N 2 ). One then finds that 



A. 



d(N-l)—2 



_ 2d 

-y 



(36) 



The d — 1 directions in which the growth factor is of or- 
der l/{na d ) contribute an amount P/(4d) to the smallest 
exponents. The growth factors due to the free flights, 
of order 2, contribute smaller amounts. If the growth 
factors had been calculated from the projection on 5v, 
the isotropic approximation would have been far less ef- 
fective, and only the growth in d — 1 directions would 
have been found. The other terms in the smallest expo- 
nents would have been absent. The same dependence on 
v would have been found, but the prefactors would have 
been smaller. 

Combining Eq. J^SJ with Eqs. and (J30J, after nu- 
merical integration over the growth factors to calculate 



10 



the averages, yields estimates for the leading order of the 
smallest exponent, at low densities, 

_ Jo.26P if d = 2 , 
Ad(JV - 1} - 2 ^ \0.32P if d = 3. (3?) 

These results are similar to the lowest exponent found 
from simulations of hard disks and hard spheres, which 
are also proportional to v [Eq. (jT5|) ]. The prefactors ob- 
tained from the estimation based on the iterative ap- 
proach and isotropic approximations differ from the sim- 
ulation results by less than 20 %. 

VI. CONCLUSIONS 

In this paper, the similarities between the chaotic prop- 
erties of a point particle colliding elastically with isotrop- 
ically distributed cylinders and systems of freely moving 
hard disks were investigated. Firstly, Monte Carlo sim- 
ulations of the spectrum of the isotropic billiard were 
presented. The lower three-quarters of the positive Lya- 
punov exponents were found to be similar to the expo- 
nents of hard disks as a function of the density and par- 
ticle number. The larger exponents behave differently, as 
was to be expected. 

Further, an analytical estimate of the smallest Lya- 
punov exponent of hard disks not due to localization was 
discussed. In this calculation, the collective property of 
the eigenvector belonging to the smallest exponent was 
used, by approximating the distribution of scatterer ori- 



entations as isotropic. This makes it possible to use the 
techniques developed in Ref. to calculate Lyapunov 
exponents through the partial stretching factor. Based 
on the calculations presented in Ref. |]|, an approxima- 
tion for the partial stretching factor was made. The re- 
sults of this calculation resemble results of the simulation 
for hard disks. The smallest exponents depend on the 
collision frequency v in the correct way and are indepen- 
dent of the particle number N, for sufficiently large N. 
The linear dependence on v of the smallest exponents is 
entirely due to the shape of the scatterers. The prefactor 
deviates from the simulation results by about 20%. 

Both calculations discussed in this paper indicate that 
the lower end of the spectrum of hard disks or spheres 
is predominantly determined by the shape of the scat- 
terers, and not so much by the scatterer orientations. 
With better approximations of the partial stretching fac- 
tor, it should be feasible to use the method developed in 
Ref. 01 with an isotropic distribution of scatterer ori- 
entations to calculate a large portion of the lower end 
of the Lyapunov spectrum of hard disks and spheres. A 
large portion of the Lyapunov spectrum of many-particle 
systems can be understood by using this isotropic ap- 
proximation. 
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